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ABSTRACT 



It has been suggested that the recently discovered exoplanet GJ581d might be able to support 
liquid water due to its relatively low mass and orbital distance. However, GJ581d receives 35% 
less stellar energy than Mars and is probably locked in tidal resonance, with extremely low 
insolation at the poles and possibly a permanent night side. Under such conditions, it is unknown 
whether any habitable climate on the planet would be able to withstand global glaciation and / 
or atmospheric collapse. Here we present three-dimensional climate simulations that demonstrate 
GJ581d will have a stable atmosphere and surface liquid water for a wide range of plausible cases, 
making it the first confirmed super-Earth (exoplanet of 2-10 Earth masses) in the habitable zone. 
We find that atmospheres with over 10 bar C0 2 and varying amounts of background gas (e.g., 
N2) yield global mean temperatures above °C for both land and ocean-covered surfaces. Based 
on the emitted IR radiation calculated by the model, we propose observational tests that will 
allow these cases to be distinguished from other possible scenarios in the future. 



Subject headings: astrobiology — planet-star interactions — planets and satellites: 
techniques: spectroscopic 



atmospheres- 



1. Introduction 

The local red dwarf Gliese 581 (20.3 ly from 
the Sun, M = 0.31Ms„ n , L = 0.0135Ls U n, spec- 
tral type M3V) ( |Hawley et~aT]|1997[ ) has received 
intense interest over the last decade due to the low 
mass exoplanets discovered around it. As of early 
2011 it has been reported to host up to six planets 



Udry et al. 


2007 


Mayor et al. 


2009 


Vogt et al. 



2010). One of these, GJ581g, was announced in 
September 2010 and estimated to be in the hab- 
itable zone (the orbital range in which a planet's 
atmosphere can warm the surface sufficiently to al- 



low surface liquid water) ( |Kasting et al. 1993 Pier- 
rehumbert||2011 ). However, its discovery has been 
strongly disputed by other researchers, including 



the team responsible for finding the other four 



planets in the system (Kerr 2010 Tuomi 2011) 



For the moment, therefore, GJ581g remains un- 
confirmed. 

GJ581d, in constrast, which was first discov- 
ered in 2007 and has a minimum mass between 
5.6 and 7.1 M Earth, has now been robustly con- 



firmed by radial velocity (RV) observations ( Udry 



et al. 2007 Mayor et al. 2009| Vogt et al. 2010 1 



Due to its greater distance from the host star, 
GJ581d was initially regarded as unlikely to have 
surface liquid water unless strong warming mech- 



anisms due to e 


g., CO2 clouds (Forget and Pier- 


rehumbert|1997 


Selsis et al.|2007 


) were present in 



its atmosphere. Recently, simple one-dimensional 



radiative-convective studies ( 


Wordsworth et al. 


2010b[ |von Paris et al.||2010[ 


Kaltenegger et al. 



2011) have suggested that a dense atmosphere 
could provide a significant greenhouse effect on 
GJ581d. However, the planet's tidal evolution 
poses a key problem for its habitability. 

As it is most likely either in a pseudo-synchronous 
state with a rotation period that is a function of 
the eccentricity, or in spin-orbit resonance like 
Mercury in our Solar System ( Leconte et al.|2 010 



Heller et al.|2011 ), GJ581d should have extremely 
low insolation at its poles and possibly a perma- 
nent night side. Regions of low or zero insolation 
on a planet can act as cold traps where volatiles 
such as H2O and CO2 freeze out on the surface. 



A few previous studies (Joshi et al. 1997 Joshi 



2003) have examined atmospheric collapse in 3D 
with simplified radiative transfer, but only for 
Earth- like atmospheric pressures or lower (0.1 to 
1.5 bar). For low values of stellar insolation and 
large planetary radii, even dense CO2 atmospheres 
will be prone to collapse, which could rule out a 
stable water cycle altogether for a super-Earth 
like GJ581d. To conclusively evaluate whether 
GJ581d is in the habitable zone, therefore, three- 
dimensional simulations using accurate radiative 
transfer are necessary. 

Here we present global climate model (GCM) 
simulations we performed to assess this issue. In 
Section [2] we describe the model used. In Section 
[3] we describe our results, and in Section [4] we dis- 
cuss implications and propose future observational 
tests for the simulated habitable scenarios. 

2. Method 

In our simulations we made the initial hypoth- 
esis that GJ581d has a climate dominated by the 
greenhouse effects of CO2 and / or H 2 0, as is 
the case for all rocky planets with atmospheres 
in the Solar System (Venus, Earth and Mars). To 
assess the influence of water on the climate in- 
dependently, we considered two classes of initial 
condition: a rocky planet with no water, and an 
ocean planet, where the surface is treated as an 
infinite water source. C0 2 was taken as the pri- 
mary constituent of the atmosphere and H2O was 
allowed to vary freely, with surface ice / liquid 
and cloud formation (including radiative effects) 
taken into account for either gas when necessary. 



Restricting the composition of the atmosphere to 
two species in this way allows us to determine con- 
servative conditions for habitability, as it neglects 
the warming due to other greenhouse gases like 
CH4 or buffer gases like 



or Ar ( 


von Paris et al. 


Li et al.|2009 


1. 



The simulations were performed using a new 
type of GCM that we developed specifically for ex- 
oplanct and paleoclimate studies. It uses radiative 
transfer data generated directly from high resolu- 
tion spectra, which allows the accurate simulation 
of climates for essentially any atmospheric cock- 
tail of gases, aerosols and clouds for which optical 
data exists. The dynamical core of the model was 
adapted from the LMD Mars GCM, which uses 
an enstrophy and total angular momentum con- 
serving finite difference scheme ( |Sadourriy 1975 



Forget et al.p 999). Scale-selective hyperdiffusion 
was used in the horizontal plane for stability. The 
planetary boundary layer was parameterised us- 



ing the method of Mellor and Yamada ( 1982 1 to 



calculate turbulent mixing, with the latent heat of 
H2O also taken into account in the surface tem- 
perature calculations when necessary. A standard 
roughness coefficient of zq = 1 x 10~ 2 m was used 
for both rocky and ocean surfaces for simplicity, 
although we verified that our results were insensi- 
tive to variations in this parameter. Spatial reso- 
lution of 32 x 32 x 20 in longitude, latitude and al- 
titude was used for all simulations; we performed 
one comparison test at the highest rotation rate 
with 64 x 64 x 20 resolution and found that the 
differences were small. 

Our radiative transfer scheme was similar to 
that we developed previously for one-dimensional 
simulations ( |Wordsworth et ah 2010a|b ). For a 
given mixture of atmospheric gases, we computed 
high resolution spectra over a range of tempera- 
tures, pressures and gas mixing ratios. For this 
study we used a 6 x 9 x 7 temperature, pressure 
and H2 O volume mixing ratio grid with values T = 
{100,150, ...,350} K, p = {10- 3 ,10- 2 ,...,10 5 } 
mbar and qH 2 = {10 -7 , 10~ 6 , . . . , 10 -1 }, respec- 
tively. The correlated-A; method was used to pro- 
duce a smaller database of coefficients suitable 
for fast calculation in a GCM. The model used 
38 spectral bands in the longwave and 36 in the 
shortwave, and sixteen points for the <?-space in- 
tegration, where g is the cumulated distribution 
function of the absorption data for each band. 



2 



In most simulations CO2 was assumed to be the 
main constituent of the atmosphere, except the 
locally habitable ice planet experiments (see Sec- 
tion[3j), where N 2 was used. CO2 collision-induced 
absorption was included using a parameterisation 
based on the most recent theoretical and experi- 



Manabe and Wetherald (1967). Precipitation of 



mental studies ( 


Wordsworth et al.| 


2010a 


Gruszka 


and Borysow 


1998 


Baranov et al. 


2004) 


. For the 



stellar spectrum, we used the Virtual Planet Lab- 
oratory AD Leo data (Segura et al. 2005). AD 



Leo is a M-class star with effective temperature 
T e ff = 3400 K, which is acceptably close to the 
most recent estimates of T e ff = 3498 ± 56 K 



von Braun et al. 2011) for the purposes of cli- 



mate modelling. A two-stream scheme (Toon et al 



1989) was used to account for the radiative ef- 



fects of both clouds and Raylcigh scattering, as in 



Wordsworth et al. (2010b) 



In the water cycle and cloud modelling, care 
was taken to ensure that the parameterisations 
used were based on physical principles and not 
tuned to Earth-specific conditions. When this 
was not possible (as for e.g., the density of con- 
densable cloud nuclei in the atmosphere N c ), we 
tested the sensitivity of our results to variations in 
those parameters. Three tracer species were used 
in our simulations: CO2 ice, H2O ice and H 2 
vapour. Tracers were freely advected in the at- 
mosphere, subject to changes due to sublimation 
/ evaporation and condensation and interaction 
with the surface. For both gases, condensation was 
assumed to occur when the atmospheric temper- 
ature dropped below the saturation temperature. 
Local mean CO2 and H2O cloud particle sizes were 
determined from the amount of condensed mate- 
rial and the density of condensable nuclei N c . This 
parameter was set to 10 5 kg -1 in most of our sim- 
ulations; we tested the effect of varying it over the 
range 10 4 to 10 6 kg -1 and found that the maxi- 
mum difference in mean surface temperature after 
60 orbits was less than 5 K. As a further test of 
the robustness of our results, we also performed 
some tests with cloud radiative effects removed al- 
together (see Section 3). 

Ice particles of both species were sedimented 
according to Stokes law ( Forget et al.|[l9 99). Be- 
low the stratosphere, adjustment was used to relax 
temperatures due to convection and / or conden- 
sation of CO2 and H2O. For H2O, moist and large- 
scale convection were taken into account following 



H 2 due to coagulation was also included using a 



simple threshold parameterisation (Emanuel and 
Ivkovi-Rothman||1999 1 . 



On the surface, the local albedo varied accord- 
ing to the composition (rocky, ocean, CO2 or H 2 
ice; see Table [IJ. In the wet simulations, ice for- 
mation (melting) was assumed to occur when the 
surface temperature was lower (higher) than 273 
K, and temperature changes due to the latent heat 
of fusion were taken into account. In all cases, 
the simulations were run until collapse/glaciation 
occurred or steady states of thermal equilibrium 
were reached. The time taken to reach ther- 
mal equilibrium can be estimated from the atmo- 



spheric radiative relaxation timescale (Goody and 
YTjngfl989 l 

CpP s 



ogTl ' 



(1) 



where c p , p s , g, T e and a are the specific heat ca- 
pacity of the atmosphere, the mean surface pres- 
sure, the surface gravity, the atmospheric emis- 
sion temperature and the Stefan-Boltzmann con- 
stant, respectively. Taking T e = 200 K, c p = 850 
J K _1 kg -1 , p s = 30 bar and g — 16.6 m s~ 2 
yields r r ~ 4000 Earth days (60 GJ581d orbits). 
This was close to the timescales we observed in 
the model by plotting time series of mean surface 
temperature. 

We used the minimum mass for GJ581d given 



by Mayor et al. ( 2009 ) instead of the smaller value 



(Mmin = 5.6M Earth ) proposed by |Vogt et al 
( |2010[ ). We took the actual mass of GJ581d to 
be M = M mm /sm 60° = 8.2M Earth , given that 
the statistically most probable value for the in- 
clination angle is 60°. Radius and gravity for 
rocky and ocean / ice cases (Table [T]) were then 
determined from theoretical models ( |Sotin et al 
2007). In the latter case, the assumed bulk com- 



position of the planet was 50% H 2 0. Model tests 
using M = M min and M = 1.6M min (the latter 
value comes from dynamical stability considera- 



tions (Mayor et al. 2009 )) did not reveal significant 
differences from our main results. 

To produce emission spectra, we recorded top- 
of-atmosphere longitude-latitude maps of outgoing 
fluxes computed by the GCM over one orbit. In- 
clination angle of the orbit relative to the observer 
was assumed to be 60° . An isotropic (Lambertian) 
distribution of specific intensities at the top of 



3 



the atmosphere was assumed; comparison with a 
line-by-line radiative transfer code at wavelengths 
where the limb-darkening was most pronounced 
revealed that the disk-integrated flux error due to 
this effect was below 5%. 

3. Results 

We performed simulations with 5, fO, 20 and 
30 bar atmospheric pressure and 1:1, 1:2 and 1:10 
orbit-rotation resonances for both rocky and ocean 
planets (see Table 1). Eccentricity was set to zero 
in the simulations we present here, with the stellar 
flux set to the true value at a = 0.22 and not in- 
creased to account for orbital averaging as in our 
ID study ( Wordsworth et al.|2010b ). However, we 
also performed tests with e = 0.38 and found that 
eccentricity was not critical to the results. Simi- 
larly, we assumed zero obliquity to assess the like- 
lihood of atmospheric collapse in the most severe 
cases, but sensistivity tests showed that its influ- 
ence on the climate was second order at the high 
pressures where the atmosphere remained stable. 

CO2 gas had a powerful greenhouse warming 
effect in our simulations because GJ581d orbits a 
red dwarf star. In the Solar System, Rayleigh scat- 
tering reflects the bluer incoming starlight much 
more effectively, and hence a planet receiving the 
same flux as GJ581d would be uninhabitable for 
any CO2 pressure. In our rocky simulations, for 
pressures below ~10 bar the atmosphere was in- 
deed unstable, and began to condense on the dark 
side and / or poles of the planet (Fig. 2a). How- 
ever, for denser atmospheres, we found that hori- 
zontal heat transport and greenhouse warming be- 
came effective enough to remove the threat of col- 
lapse and allow surface temperatures above the 
melting point of water (Fig. 1). Rotation rate 
affected the atmospheric stability through both 
the insolation (the planet has permanent day and 
night sides in the most extreme 1:1 resonance case) 
and the horizontal heat transport (faster rotat- 
ing atmospheres were less efficient at transporting 
heat polewards). 

Even in the stable simulations, some CO2 usu- 
ally condensed in the middle atmosphere (Fig. 3) , 
leading to CO2 ice cloud formation as occurs on 



sumptions, due to the scattering of infrared radi- 



ation from the ground ( Forget and Pierrehumbert 



1997). To check the robustness of our results, we 



also performed tests with the (unrealistic) assump- 
tion of no CO2 cloud radiative effects, and found 
that the atmospheres were colder, but still stable 
above around 10 bar. 

In the ocean planet simulations, we neglected 
all oceanic horizontal heat transport. This as- 
sumption led us to overestimate global tempera- 
ture differences and hence the probability of atmo- 
spheric collapse, in keeping with the aim of a con- 
servative habitability estimate. We found that at- 
mospheric H 2 vapour greatly increased warming, 
while H2O cloud formation low in the atmosphere 
tended to cool the planet by increasing the plan- 
etary albedo (Fig. 3). This led to a transition in 
the climate as total pressure increased. At 10 bar 
and below, cooling effects dominated and runaway 
glaciation occurred, followed by atmospheric col- 
lapse. From 20 bar, the positive feedback of water 
vapour on greenhouse warming raised mean tem- 
peratures significantly compared with the rocky 
case (Fig. 2b) . In the intermediate region, climate 
stability was difficult to assess as surface temper- 
ature trends were extremely small (as low as a few 
K per 100 orbits in some cases). There, our re- 
sults were somewhat sensitive to our microphysical 
assumptions; for example, lowering the precipita- 
tion threshold Iq (see Table 1) resulted in optically 
thinner H2O clouds and stable, cooler climates 
at lower total pressures. Nonetheless, at higher 
pressures we found stable, hot climates with little 
global temperature variations even in the tidally 
locked cases, regardless of the choice of microphys- 
ical parameters. 

We investigated the effect of starting with a 
surface covered by H2O ice. For the dense atmo- 
spheres where the climate was stable, this only 
weakly changed the planetary albedo, and the ice 
melted for pressures of 20 bar or more. Hence 
the runaway (H2O) glaciation that is possible on 
Earth (Budyko 1969[ ) would be unlikely to occur 
on an ocean-planet GJ581d. Our results here con- 



present-day Mars (Montmcssin et al. 2007). CO2 



clouds typically had a net warming effect (up to 
12 K), depending on the cloud microphysical as- 



trast with those of Spiegel et al. ( 2009[ ), who used 
a ID energy balance model to investigate the ef- 
fects of obliquity on planetary habitability, be- 
cause they neglected the dependence of planetary 
albedo on the total pressure in dense atmospheres. 
We also investigated the possibility that an ice- 



4 



covered, tidally locked GJ581d with permanent 
day and night sides could be locally habitable on 
the day side due to partial melting. However in 
this scenario, the day side only stayed warm when 
the atmosphere was thin and hence inefficient at 
transporting heat. As a result, the planet's dark 
side became cold enough for the collapse of even 
an N 2 atmosphere. An ice planet with only a thin 
sublimation-driven H 2 atmosphere could have 
dayside temperatures above 273 K and contin- 
ual transport of H2O to the dark side, but the 
low atmospheric pressure would preclude liquid 
water except in extremely limited sub-surface re- 
gions. Clearly, dense, stable atmospheres offer 
better prospects for habitability. 

4. Discussion 

Are CO2 partial pressures of over 10 bar a real- 
istic possibility for GJ581d? Given the planet's as- 
sumed gravity, this corresponds to a CO2 column 
of only 4 to 6 bars on Earth. Even though Venus 
and Earth are smaller, their total CO2 invento- 
ries scaled by planetary mass are around 10-100 
times greater than this. On Venus, it is thought 
that a large fraction of this inventory is in the at- 
mosphere ( Bullock and Grinspoon|1996 ) , while on 
Earth, the atmospheric partial pressure is regu- 
lated over geological timescales by the carbonate- 
silicate cycle ( Walker et al.|1981 1 via plate tecton- 
ics. The geophysics of GJ581d is unknown, but if 
a similar mechanism were present there, its atmo- 
spheric CO2 would stabilize above the level needed 
to maintain a liquid water cycle by negative feed- 
back. 

C02-H 2 0-(buffer gas) atmospheres are clearly 
only a subset of all the possible scenarios for 
GJ581d. If small quantities of additional green- 
house gases such as CH4 or SO2 were present, 
they would increase warming further. However, if 
these gases were chemically unstable, they would 
need to be continually emitted by e.g., volcanic 
or biological sources to have a long-term effect on 
the climate. More drastically, the planet could 
have a thick H2-He envelope like Uranus or Nep- 
tune, or its atmosphere could have been removed 
entirely by increased stellar activity early in the 
life of the host star. Both scenarios would ex- 
clude surface liquid water (the former through an 
excess of surface pressure). GJ581d's large mass 



means that it can retain an atmosphere more eas- 
ily than Earth or Venus, but the increased XUV 
and ion fluxes from M-class stars early in their his- 
tory means that removal processes can be much 



more efficient than in the Solar System ( Tanigawa 



and Ikoma 


2007 


Tian 


2009) 



The escape of at- 
mospheric hydrogen depends sensitively on both 
XUV flux levels and radiative cooling by Hjj" ions, 
both of which are poorly constrained for GJ581d 
dSelsis et al.||2007l iKoskinen et al.||2007b. More- 



over, stripping by the stellar wind, which does not 
discriminate between light and heavy atoms, could 
have removed large amounts of any gas (including 
C0 2 ) (Lammer et al. 20071. Detailed modelling 
of atmospheric escape may be able to constrain 
the effects of these processes on the atmosphere 
further. 

Unlike the majority of the Kepler planetary 
candidates, Gliese 581d is relatively close to Earth, 
so in the future it will be possible to establish 
which scenario applies to it through direct spectro- 
scopic observations. To determine the sensitivity 
required for this, we used the outgoing longwave 
radiation (OLR) from our simulations to produce 
synthetic emission spectra for the habitable cli- 
mates just discussed, along with those for a planet 
with no atmosphere (Fig. 4) . We did not calculate 
emission spectra for the H 2 -He dominated case, 
but equilibrium chemistry calculation^] allowed us 
to establish that for an atmosphere where H 2 is 
the dominant species by mass, the abundance of 
CO2 should be below the limit of dctectability. 
For a rocky / icy planet where volatiles like C0 2 
have collapsed on the surface and the atmosphere 
is thin or non-existent, the flux variations over one 
orbit will be large (grey region in Fig. [4]). Obser- 
vations of CO2 and H2O absorption bands with 
low phase variations (red / blue regions in Fig. 
[4]) would hence be a strong indicator of the kind 
of stable, dense habitable atmospheres discussed 
in this paper. From Fig. [4j the flux sensitivity 
necessary for this in the infrared will be of order 
10~ 21 W raT 2 /im _1 , corresponding to a planet 
/ star ratio lower than 10~ 6 for wavelengths < 
20 y«m. This is beyond the capabilities of current 
space and ground-based observatories, but it will 



1 These were performed using a code developed by R. 
Bounaceur (LRGP, CNRS, Nancy, Prance), based on an 
algorith m to minimize t he free energy of a given mixture 
of gases | |Reynolds|1986| . 
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become possible with future improvements in in- 
strumental technology. 

F.S. acknowledges support from the Euro- 
pean Research Council (Starting Grant 209622: 
E 3 ARTHs). 
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Table 1: Standard parameters used in the climate 
simulations. 



Stellar luminosity 
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r t l 


0.U135 


Orbital semi-major axis 
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Orbital eccentricity 
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Obliquity 
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Initial atmospheric pressure 


Ps 


[bars] 


5-30 


Radius (rocky) 


r 


[^Earth] 


1.8 


Radius (ocean) 


r 


l^Earth] 


2.3 


Surface gravity (rocky) 
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[m S~ 2 ] 


25.0 


Surface gravity (ocean) 
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[m S -2 ] 


16.6 


Surface albedo (rocky) 


A, 




0.2 


Surface albedo (ocean) 


A ; 




0.07 


Surface albedo (ice) 


A 
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Surface roughness coefficient 




\m] 


1 x 10- 


Precipitation threshold 




[kg kg *] 


0.001 


Number of cloud condensation nuclei 


N, 


3 [kg" 1 ] 


1 x 10 5 



This 2-column preprint was prepared with the AAS IATgX 
macros v5.2. 
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Fig. 1. — Surface temperature snapshots after 60 
orbits integration time for rocky planet simula- 
tions with a) 1:1, b) 1:2 and c) 1:10 tidal resonance 
and a 20-bar CO2 atmosphere. 
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Fig. 3. — Latitude-altitude map of H2O (blue) and 
CO2 (red) cloud coverage for a 20-bar ocean sim- 
ulation with 1:2 resonance, averaged in longitude 
and over one orbit. While the cloud deck altitudes 
were similar for all simulations, the latitudinal dis- 
tribution depended strongly on the atmospheric 
dynamics and hence on the rotation rate and to- 
tal pressure. 
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Fig. 2. — Surface temperature vs. pressure for 
a) rocky and b) ocean planet simulations. Red, 
green and blue markers represent tidal resonances 
of 1:1, 1:2 and 1:10. Points are separated for 
clarity; simulations were performed at 5, 10, 20 
and 30 bar for each all cases. Crosses and error 
bars show mean and maximum / minimum tem- 
peratures (sampled over one orbit and across the 
planet's surface) in stable cases, while circles in- 
dicate simulations where the atmosphere began to 
collapse or runaway glaciation occurred. 
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Fig. 4. — Infrared emission spectra for an ocean 
/ rocky GJ581d with a 20-bar atmosphere (blue 
/ red) and a rocky GJ581d with no atmosphere 
(grey). Thickness of the lines corresponds to the 
maximum / minimum difference in flux over one 
orbit, for the most probable observation angle of 
60°. The maximum (minimum) surface temper- 
ature in the airless case is 271 K (37 K). For 
the cases with atmospheres, emission of the ocean 
planet is higher at most wavelengths due to the 
increased planetary radius. The major absorption 
features are labelled by molecule / process (C0 2 - 
CO2 corresponds to CO2 collision-induced absorp- 
tion). 
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